home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / findiff.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.1 KB  |  43 lines  |  [MATF/MATL]
  1. function [T,X] = findiff(p,q,r,a,b,alpha,beta,n)
  2. % [T,X] = findiff(p,q,r,a,b,alpha,beta,n)
  3. % Finite difference solution to the boundary value problem
  4. % x`` = p(t)x`(t)+q(t)x(t)+r(t), x(a) = alpha, x(b) = beta
  5. % This subroutine uses trisys.m to solve a tri-diagonal system.
  6. % p  is a function, input.
  7. % q  is a function, input.
  8. % r  is a function, input.
  9. % a  is the left endpoint, input.
  10. % b  is the right endpoint, input.
  11. % alpha is the left boundary value, input.
  12. % beta is the right boundary value, input.
  13. % n  is the number of steps, input.
  14. % T  is the vector of abscissas, output.
  15. % X  is the vector of ordinates, output.
  16. T  = zeros(1,n+1);
  17. X  = zeros(1,n-1);
  18. Va = zeros(1,n-2);
  19. Vb = zeros(1,n-1);
  20. Vc = zeros(1,n-2);
  21. Vd = zeros(1,n-1);
  22. h = (b - a)/n;
  23. for j=1:n-1,
  24.   Vt(j) = a + h*j;
  25. end
  26. for j=1:n-1,
  27.   Vb(j) = -h^2*feval(r,Vt(j));
  28. end
  29. Vb(1)   = Vb(1)   + (1 + h/2*feval(p,Vt(1)))*alpha;
  30. Vb(n-1) = Vb(n-1) + (1 - h/2*feval(p,Vt(n-1)))*beta;
  31. for j=1:n-1,
  32.   Vd(j) = 2 + h^2*feval(q,Vt(j));
  33. end
  34. for j=1:n-2,
  35.   Va(j) = -1 - h/2*feval(p,Vt(j+1));
  36. end
  37. for j=1:n-2,
  38.   Vc(j) = -1 + h/2*feval(p,Vt(j));
  39. end
  40. X = trisys(Va,Vd,Vc,Vb);
  41. T = [a,Vt,b];
  42. X = [alpha,X,beta];
  43.